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ABSTRACT 

We consider the three-dimensional instability of a layer of horizontal magnetic field in 
a polytropic atmosphere where, contrary to previous studies, the field lines in the initial 
state are not unidirectional. We show that if the twist is initially concentrated inside the 
unstable layer, the mod ifications of the instability reported by several authors (see e.g. 
Cattaneo et al.| \ 199 0a)) are only observed when the calculation is restricted to two di- 
mensions. In three dimensions, the usual interchange instability occurs, in the direction 
fixed by the field lines at the interface between the layer and the field-free region. We 
therefore introduce a new configuration: the instability now develops in a weakly mag- 
netised atmosphere where the direction of the field can vary with respect to the direction 
of the strong unstable field below, the twist being now concentrated at the upper inter- 
face. Both linear stability analysis and non-linear direct numerical simulations are used 
to study this configuration. We show that from the small-scale interchange instability, 
large-scale twisted coherent magnetic structures are spontaneously formed, with possi- 
ble implications to the formation of active regions from a deep-seated solar magnetic 
field. 
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1 Introduction 

Active regions at the surface of the Sun are believed to be the vis- 
ible manifestation of deep-seated intense magnetic fields. The 
generation of these predominantly toroidal fields is associated 
with the differential rotation in the tachocline. However, the 
transport of these strong fields from the lower part of the con- 
vective zone up to the photosphere remains one of the major 
unknowns of the solar cycle. It is clear that super-equipartition 
fields are buoyant and could therefore rise up to the surface, but 
the observed size, coherence and twist of the magnetic structures 
forming active regions remain largely unexplained. 

In order to obtain some understanding of the transport of 
magnetic structures by buoyancy, highly idealised models have 
been considered. An early model ([Parker 1955) considers the 
buoyancy properties of the magnetic field when it is concen- 
trated in slender pressure-confined flux tubes. A simple model 
describing the dynamical evolution of such structures was pro- 
posed by Spruit (198 lj and was subsequently used by several 
authors (see for example Moreno-Insertis 1983). This model has 
successfully reproduced Joy's law for the tilt of bipolar regions 
at the solar photosphere (Choudhuri 1987 ; D'Silva & Choudhuri 
|1993| [Caligari et al.|1995) as well as the relationship between the 
tilt angle and the total flux of active regions ( [Fan et al.|1994| ). 
Many authors have sought to move on from this simple model 



Email address: b.favier@damtp.cam.ac.uk 



and have considered the full set of compressible magnetohydro- 
dynamic (MHD) equations (or the anelastic approximation) to 
study the evolution of a simplified magnetic structure mimick- 
ing a flux tube. These models have shown that a certain amount 
of twist must be present inside the tube in order for it to remain 
coherent as it rises; such a twist corresponds to a non-zero axial 
current inside the flux tube. Without twist, the flux tube is shred- 
ded apart by the vorticity generated by its own motion ( Schiissler 



1979 



Longcope et al.|1996]|Emonet & Moreno-Insertis [1996 



Fan et al. ||1 998 [ | Jouve & Brun 2007). Even more problematic 



than the issue of tube coherence is the interaction between the 
rising structures and the small-scale convective motions. Some 
studies have considered the rise of gradually twisted magnetic 
flux tubes through the convective zone both in Cartesian ([Fan 
|et~aLl[2003l |Abbett et alj|2004] ) and spherical ( |Jouve & Brun 



2009) geometries. Again, a given amount of twist is necessary 
for the flux tubes to remain coherent and rise through the convec- 
tive zone. There is observational evidence of non-potential mag- 
netic fields in active regions (Seehafer 1990, Leka et al. 1996; 
|Pevtsov et al.|2Q0l) , with the implication that the flux tubes are 
twisted before they emerge as active regions. While the twist is 
no more than one turn in general (Chae & Moon 2005), the cor- 
responding Lorentz forces are implicated in chromospheric flux 
eruptions (Torok & Kliem 2005). The origin of the twist is still 
unclear as there are various different explanations. The effect 
of the Coriolis force (Fan & Gong 2000) and differential rota- 
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tion ([DeVore 2000 ) on the rising structures can only contribute a 
small fraction of the observed twist values ( Holder et al. 2004). 

Other mechanisms have been suggested to explain the ori- 
gin of non-potential magnetic fields. One possibility is that the 
tube becomes twisted due to the effect of small-scale helical mo- 
tions in the convective zone acting on the rising magnetic field, 
the so-called E-effect ( [Longcope et al.|19 98). Another sugges- 
tion, which is closely related to the content of this paper, is that 
the twist arises by accretion of poloidal fields during the rise of 
a flux tube ( |Choudhuri|2003||Chatterjee et al.|2006] ). However 
indirect observational evidence suggests the presence of an ini- 
tial twist in flux tubes at the base of the solar convective zone 
( [Holder et al.|2004] >, so that the twist may be produced during 
the formation process of the flux tube itself. 

Another issue, which we don't yet fully understand, is the 
mere existence of localised magnetic structures similar to flux 
tubes below the convective zone. What is known is that the dif- 
ferential rotation profile produces a strong toroidal field at the 
base of the solar convective zone. We can attempt to model 
the appearance of isolated flux structures by looking at a uni- 
form layer of toroidal field where the vertical structure of the 
magnetic layer is initially imposed; then the horizontal length 
scale will be set naturally by the buoyancy instability. As a re- 
sult of the magnetic pressure, the density of polytropic atmo- 
sphere in the magnetised region is reduced, assuming the tem- 
perature is not modified by the presence of the magnetic field. 
Consequently, the upper interface of the magnetic slab is sus- 
ceptible to Rayleigh-Tay lor- type instabilities. This paper is prin- 
cipally concerned with the three-dimensional evolution of these 
instabilities when the initial magnetic field has non-zero twist. 
The two-dimensional evolution of states with no initial twist 
was investigated by Cattaneo & Hughes ( 1988J; since no vari- 
ation was permitted along the original field direction this treat- 
ment could only encompass the interchange instability, with the 
field lines remaining straight. This configuration has also been 
considered numerically in three dimensions by several authors 
(Matthews et al. 1995 ; Wissink et al. 2000), whereas alternative 
initial conditions (in which the temperature profile is altered by 
the magnetic field instead of the density) were considered by 
|Fan| pOOl). It has proved very difficult to produce flux struc- 
tures of significant size from this type of instability. Whatever 
the scale of the original instability secondary vortex instabilities 
due to the down flows between the plumes disrupt the magnetic 
structures leaving a rather diffuse and dynamically inactive field 
. Note that an alternative model has been suggested by |KersalI| 
|et al.| ( [2007] ), where they considered the non-linear evolution of 
magnetic buoyancy instabilities resulting from a smoothly strati- 
fied horizontal magnetic field. Cattaneo et al. (1990a) conducted 
two-dimensional simulations of an initial twisted field structure. 
They found that in their constrained geometry large structures 
could appear as result of the twist. We will show in this paper 
that this effect disappears if we allow the motion to be fully 
three-dimensional. 

The earlier simulations to examine the evolution of buoy- 
ant magnetic structures assumed that the region above the initial 
magnetic slab was field free. In this case the initial potential en- 
ergy contained in the initial condition is eventually transferred to 
kinetic energy and is finally dissipated. Any possible effects of 
the magnetic field existing above the unstable slab are neglected. 
Is this a reasonable simplification? 

As shown by many simulations of overshooting convection 
in the presence of magnetic field (see for example Tobia s et ah] 
( |2001| )), the strong toroidal field initially injected in the system 
or generated by a shear flow mimicking the tachocline ( Guerrero 
|& Kapyla|2011| ) is predominantly contained in the stable radia- 
tive zone. However, despite the effects of turbulent pumping, a 



non-negligible part of the field is still present in the convection 
zone, both due to the redistribution of large-scale magnetic fields 
and to local small-scale dynamo action. If the strong toroidal 
field is buoyantly unstable, it will rise and interact with the small 
scale magnetic perturbations present in the convective zone. 

The purpose of this paper is to investigate the simplest 
model that allows us to consider the buoyancy instability of 
a toroidal magnetic field in a weakly magnetised atmosphere 
above. The layer of strong toroidal field models the field induced 
by shear in the tachocline whereas the magnetised atmosphere 
above represents the magnetic perturbations in the convective 
zone. For this first, illustrative study, the field in the convective 
zone is assumed to be unidirectional but with a different direc- 
tion as in the layer of strong field below. This is of course highly 
idealised, but allows us to derive a model involving few free pa- 
rameters. The instability of a sheared magnetic field has already 
been studied ( [Cattaneo et al.| ( |1990a1b| ); |Kusano et al.| ( |1998| »; 
|Nozaw a ( 2005 )). In this case, the twist is uniformly distributed 
inside the unstable slab and the atmosphere above is field-free. 
To our knowledge, this type of instability has only been consid- 
ered in published papers using two-dimensional numerical sim- 
ulations. The configuration we consider later in the paper is dif- 
ferent as the field lines are unidirectional except in a thin region 
where the twist and current are concentrated. A similar setup was 
considered by |Stone & Gardiner ( 2007b a) where they look at 
the effect of a transverse magnetic field on the Rayleigh-Taylor 
instability. 

The paper will proceed as follows: in section|2] we describe 
the governing equations, the model and physical parameters and 
the numerical methods. We then extend the results by [Cattaneo] 
|et al.| ( | 1990a] ) and |Kusano et al.| ( [T998] ) to the three-dimensional 
case in section [3] Finally, the interchange instability in a weakly 
magnetised atmosphere is considered in section [4] using both 
linear stability analysis and non-linear numerical simulations in 
three dimensions. 



2 Model and method 

2.1 Model and governing equations 

We consider the evolution of a plane-parallel layer of compress- 
ible fluid, bounded above and below by two impenetrable, stress- 
free boundaries, a distance d apart. The upper boundary is held 
at fixed temperature, T whereas a vertical temperature gradient 
9 is fixed at the lower boundary. The geometry of this layer is 
defined by a Cartesian grid, with x and y corresponding to the 
horizontal coordinates. The z-axis points vertically downwards, 
parallel to the constant gravitational acceleration g = gz. The 
horizontal size of the fluid domain is defined by the aspect ratios 
X x and X y , so that the fluid occupies the domain < z < d, 
< x < X x d and < y < X y d. The physical properties of 
the fluid, namely the specific heats c p and c v , the shear viscos- 
ity fi, the thermal conductivity K, the magnetic permeability /io 
and the magnetic diffusivity 77, are assumed to be constant. The 
model is identical to that used by, for example, Matthews et al. 
( fT^ , |Sl!vel^^ and |Favier & Bushby| ( |2012| )~ 

It is convenient to introduce dimensionless variables, so 
we adopt the scalings described in Matt hews et al.| ( |1995| > and 
Bushby et al. (2008). Lengths are scaled with the depth of the 
layer, d. The temperature, T, and the density, p, are scaled with 
their values at the upper surface, To and po respectively. The 
velocity, u, is scaled with the isothermal sound speed, y/R^To, 
at the top of the layer, where i?* is the gas constant. We adopt 
the same scaling for the Alfven speed, which implies that the 
magnetic field, B, is scaled with \/popoR^T . Finally, we scale 
time by an acoustic time scale dj \/R*Tq. 
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We now express the governing equations in terms of these 
dimensionless variables. The equation for conservation of mass 
is given by 

d P 



Table 1. Parameters values considered in this paper 



dt 



+ V • (pu) = o . 



(i) 



Similarly, the dimensionless momentum equation can be written 
in the following form, 

^ + u • Vu = --VP + - (V x B) x B+ 

at p p 

0(m + l)z+— V-t, (2) 
P 



where P is the pressure given by the equation of state P 
and t is the rate of strain tensor defined by 



dxj 



+ 



duj_ 
dxi 



duk 
' dx k 



pT 



(3) 



Several non-dimensional parameters appear in equation (|2| in- 
cluding 0, the vertical temperature gradient at the lower bound- 
ary. Without magnetic field, it is initially equal to AT /To where 
AT is the temperature difference between the upper and lower 
boundaries, m = gd/R*AT — 1 corresponds to the poly- 
tropic index. The dimensionless thermal diffusivity is given by 
k = K /dpoCp(R^To) 1 ^ 2 and a = pc p /K is the Prandtl num- 
ber. The dimensionless strength of the magnetic field is defined 
by F = Bq /(popoR*To), where Bo is the amplitude of the im- 
posed magnetic field. It is related to the plasma f3 by F = 2/(3. 
The induction equation for the magnetic field is 



OB 

dt 



V x (it x B - Co^V x B) 



(4) 



where Co = r\c v po /K is the ratio of magnetic to thermal dif- 
fusivity at the top of the layer (inverse Roberts number). The 
magnetic field is solenoidal so that 



V-B = 0. 



(5) 



Finally, the heat equation is 

^ + u • VT = - ( 7 - 1) TV • u + ^ V 2 T+ 
ot p 

^L^(aT 2 /2 + (oF\VxB\ 2 ) , (6) 
P 

where 7 — c p /c v . 

To complete the specification of the model, some bound- 
ary conditions must be imposed. In the horizontal directions, all 
variables are assumed to be periodic. As has already been de- 
scribed, the upper and lower boundaries are assumed to be im- 
permeable and stress-free, which implies that u XjZ = u y , z — 
u z — at z — (the upper boundary) and z — 1 (the lower 
boundary). Having non-dimensionalised the system, the ther- 
mal boundary conditions at these surfaces correspond to fixing 
T = 1 at z = and d z T = 9 at z = 1. For the mag- 
netic field boundary conditions, we choose appropriate condi- 
tions for perfectly-conducting boundaries, which implies that 



B z — B x 



B v 



at z — and z 



2.2 Initial conditions and parameters 

In this paper, we investigate the instability of a layer of purely 
horizontal magnetic field. The vertical component of the mag- 
netic field B z is therefore initially zero and we define a ver- 
tical profile for the initial horizontal field Bn(z,t = 0) = 
y/B%(z) + B$(z). In this paper, we consider two main initial 
magnetic field configurations. In section|3] the horizontal field is 
non-zero only in a prescribed layer zt < z < 1, where z t is the 



Variable 



Parameter 



Value 



a 

Co 

e 

m 
F 

\ x x X v 



Thermal diffusivity 
Prandtl number 
Inverse Roberts number 
Thermal stratification 
Polytropic index 
Dimensionless magnetic field strength 
Horizontal aspect ratio 



5 x 10- 4 
5 x 1CT 2 
5 x 10- 2 
2 

1.6 
0.1 

2x8 



depth of the unstable interface and z = 1 is the lower bound- 
ary of the numerical domain. In section [4] we include a weak 
horizontal field above the strong horizontal field in the magnetic 
layer to give a magnetised atmosphere throughout the domain. 
This new configuration allows the magnetic tension to alter the 
development of the instability. For these prescribed initial mag- 
netic fields an equilibrium static state (it = 0) can be found 
by solving equations (2} and (6| for the density and tempera- 
ture profiles respectively. The temperature gradient is slightly 
increased from the polytropic atmosphere due to the presence of 
ohmic heating, whereas the density profile is modified in order 
to keep the interface in pressure equilibrium (as seen later in fig- 
ure [2|b)), making the layer of strong field buoyantly unstable. 
This magnetised polytropic layer, coupled with a small random 
thermal perturbation, is an appropriate initial condition for our 
numerical simulations. 

There are many dimensionless parameters in this system 
and it is not viable to complete a systematic survey of the whole 
parameter space. Throughout this paper, the polytropic index is 
fixed at m = 1.6 whereas the ratio of specific heats is given 
by 7 = 5/3 (the appropriate value for a monatomic gas). This 
ensures that the initial polytropic atmosphere is stable to con- 
vective motions (even with the increase in the temperature gra- 
dient due to ohmic heating) and any resulting instability is due 
to magnetically induced buoyancy effects. It is also an appro- 
priate value for the transition between the nearly-adiabatic solar 
convective zone and the stably- stratified radiative zone below, 
which is the main focus of the present model. 

We fix the kinetic Prandtl number to be a = 0.05 and the 
inverse Roberts number to be Co = 0.05 for the work presented 
in this paper. Note that tachocline values are much lower than 
those that can be realised numerically. These two values are 
slightly larger than the ones observed in the literature. This is 
mostly justified by numerical reasons (larger values of a and 
Co allow for larger time steps). However, we have run addi- 
tional simulations using smaller values for a and Co (down to 
a — Co = 10 -2 ), along with cases with a < Co (as it should be 
for the tachocline), and the results are qualitatively unchanged. 
Note we did not explore the possibility of the double diffusive 
instability ( Silvers et al. 2009 2010), which is much more de- 
manding numerically and will be followed up at a later point. 
The thermal stratification is fixed to be 6 — 2, which corre- 
sponds to a moderately stratified layer (the layer contains ap- 
proximately three pressure scale heights). 

As will be explained later, increasing the stratification is 
not useful in the present model. We do not vary the dimension- 
less strength of the magnetic field in this paper and so F is fixed 
to be F — 0.1. However we note here that we also tried to vary 
this parameter without any qualitative changes in the results. The 
value of F in the solar tachocline is expected to be much smaller 
than the one considered here, but this regime is challenging nu- 
merically and is therefore left for future studies. 

The numerical domain is chosen to be elongated in the y 
direction in order to accommodate for the long wavelength mod- 
ulation of the secondary instability in this direction, as reported 
by |Matthews et al.| ( |l995| >. In the transverse direction, the numer- 



© 2012 RAS, MNRAS 000,[T]{TT] 



4 B. Favier et al. 



B x (t = 0)=0 




1 .5 2.0 



X 



3D 



- 0.06 

- 0.04 

- 0.02 



B x (t = 0)^0 




Figure 1. Comparison between two and three dimensional simulations. Top line: the initial field is unidirectional, directed perpendicularly to the plane. 
Bottom line: the initial magnetic layer is sheared. On the left column, we show contours of B y from the two-dimensional simulations, whereas the 
middle column corresponds to one particular (x, z) plane inside the three-dimensional domain. On the right column, we show a volume rendering of the 
magnetic energy from the 3D simulations at the same time. All visualisations are performed at the same time t « 20, except for the bottom left figure, 
which corresponds to t « 50 (due to the fact that the growth rate of the instability is strongly reduced in this particular case see, for example, Cattaneo 
[etaTl(1990a) ). 



ical domain is large enough to accommodate approximately 20 
most unstable wavelengths in the case without transverse mag- 
netic field. A summary of our choice of parameters is given in 
Table[T] 



2.3 Numerical method 

The given set of equations is solved using a modified version of 
the mixed pseudo- spectral/finite difference code originally de- 
scribed by Matthews et al. (1995). Due to periodicity in the hor- 
izontal direction, horizontal derivatives are computed in Fourier 
space using fast Fourier transforms. In the vertical direction, a 
fourth-order finite differences scheme is used and an upwind 
stencil is applied for the advective terms. The time- stepping is 
performed by an explicit adaptive third-order Adams-Bashforth 
technique, with a variable time- step. The standard resolution is 
256 grid-points in each horizontal direction and 240 grid-points 
in the vertical direction. We also consider quasi two-dimensional 
simulations for which the variations along the y direction (along 
the field lines of the unstable layer) are neglected. A poloidal- 
toroidal decomposition is used for the magnetic field in order to 
ensure that the field remains solenoidal. 

A linear stability analysis of the equilibrium state is also 
performed in section [4?7] to determine how the initial buoyancy 
instability of the strong toroidal field B y is affected by the ad- 
dition of the transverse field B x . The system of equations |T}- 
{6} is perturbed, linearised and the ideal gas law and the mag- 
netic solenoidality condition are used to eliminate the pressure 
and one component of the fluctuating magnetic field. All per- 
turbations are of the form exp (ik x x + ik y y + st), where s is 
the complex growth rate. This gives rise to the a seventh order 
system of equations for the fluctuating variables which may be 
numerically approximated by 

SXri — A n Xn , (7) 

where x ji is the vector of discretized eigenvectors and A n is 



the matrix of discretized differential operators, using a fourth- 
order finite difference approximation. The boundary conditions 
are applied through the relevant matrix coefficients using finite 
difference schemes. The resulting matrix equation is then solved 
for the eigenvalues s and the corresponding eigenfunctions Xn 
using the Linear Algebra PACKage (LAPACK). 



3 Sheared magnetic field 

In this section, we extend the configuration already considered 
by |Cattaneo et al.|( 1990a] ) and |Kusano et al.|(1998) to the three- 
dimensional case. We consider two different cases: in the first 
one, the initial field is unidirectional, say B y , whereas in the sec- 
ond case, a horizontal transverse field B x is added. For these two 
different cases, the horizontal magnetic field Bh is the same and 
we compare both two-dimensional (2D) and three-dimensional 
(3D) simulations. By two-dimensional simulation, we mean that 
all the variables are invariant in one of the horizontal direction 
(namely in the y direction). The basic state is defined as follows: 
the initial horizontal field Bh(z, t = 0) only depends on z and 
is the same in both cases, 



B H (z) = y/Bl + Bl = \ [l.O + erf (^)] 



(8) 



where z t is the depth corresponding the top of the unstable layer 
whereas A z corresponds to the width of the transition between 
the layer and the non-magnetised atmosphere above. This choice 
is compatible with the boundary conditions defined in section [2] 
In the following, we choose zt — 0.6 and Az — 1/60. The first 
reference case is obtained by taking B x = so that the field 
is initially unidirectional, in the y direction. In the sheared case, 
and following Cattan eo" et al.| ( |1990a| ), the B x is chosen to be 
linearly dependent on the depth z according to 



B x 



r) 



if Z t <Z<Z r 

otherwise 



(9) 
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Figure 2. Initial condition for the horizontal components of the magnetic 
field, the temperature, T, and the density, p. In the top image, the dashed 
lines correspond to the vertical profile of B x and B y for the particular 
case a = 7r/4. In the bottom image, the dashed lines represent the basic 
state of temperature and density in the absence of magnetic field. The 
parameter F has been increased from F = 0.1 (the value used in this 
paper) to F = 1 in order to enhance the modifications of the basic field- 
free poly tropic state. 



(this profile is actually smoothed using error function to avoid 
spurious numerical effects due to discontinuities). The depth z r 
at which the transverse field goes to zero is called the resonant 
layer by Cattaneo et al.| |l990a| ). We choose here z r — 0.8. 
The initial amplitude of the transverse field is derived from e = 
—0.75 (the maximum value of B x initially occurs at z = 0.6 
and is B x = 0.15). Note that we tried different configurations 
with different resonant layers and different amplitudes for the x 
field, but the results are qualitatively similar. 

From this initial condition, equations (T}-(6} are solved us- 
ing the numerical method described in section [23] The simu- 
lations are run until the upper and lower boundary conditions 
have a significant effect on the generated flow. We note here 
that throughout all the results, the time is presented in terms of 
sound crossing time units. We present in figure [T] a comparison 
of two-dimensional and three-dimensional simulations for the 
two cases with and without transverse field. The cases with a 
unidirectional field B y are shown in the upper part of the fig- 
ure, whereas the cases where B x ^ initially are shown in the 
lower part. It can be seen that in the unidirectional case, both 2D 
and 3D simulations are very similar. The main difference is due 
to the lack of vortex stretching in the 2D simulation, leading to 
different mixing properties. The striking similarity between 2D 
and 3D simulations is a confirmation of the interchange nature 



of the initial instability, which is purely two-dimensional. On 
the other hand, in the sheared case (i.e. when B x ^ initially, 
lower line of figure [TJ, a drastic difference is observed between 
2D and 3D simulations. In 2D, we recover the large-scale insta- 
bility already observed by several authors (see for example [Cat- 1 
|taneo et al.| ( |1990a| ) and |Kusano eTaL] ( [T998] )). The small-scale 
interchange instability observed in the unidirectional case is al- 
tered by the tension due to the transverse field, leading to a long 
wavelength instability. This strong modification is not observed 
in 3D, where the usual interchange instability occurs. The sim- 
ilarity between the 3D simulations with and without B x field 
is clear. Figure [T] also shows a three-dimensional visualisation 
comparing the unidirectional and sheared cases using the soft- 
ware VAPOK[tl( Clyne et al.|2007| ). It is clear that the nature of 
the instability is similar in both cases, apart from the change in 
direction of the field lines. In 2D, due to geometrical constraints, 
it is not possible for the x field to undergo interchange instabil- 
ity, whereas it is possible in 3D. The results are similar for many 
different configurations (e.g. changing zt, z r ,e and the shape of 
the B x field from linear to quadratic). Thus we conclude from 
what we have observed that imposing a transverse field inside 
the unstable layer is not enough to alter interchange instability 
in three dimensions, as was previously thought using 2D simu- 
lations. We cannot rule out the possibility that this conclusion 
depends on the parameters and initial conditions considered, but 
we simply conclude that it might not be as easy as in the 2D 
case. 

In the next section, we examine another initial configura- 
tion, which leads to a significant modification of the instability, 
while working both in two and three dimensions. 



4 Instability in a magnetised atmosphere 

In this section, unlike in section [3] and previous studies, the at- 
mosphere above the layer of strong field is now also magnetised. 
Initially, we fix the vertical field to be zero everywhere whereas 
the horizontal field is only varying in the vertical direction ac- 
cording to 



B H = yjBl + Bl = 0.2 + 0.5 [l.O + erf (^^)] 



(10) 



where z t is the depth corresponding the top of the unstable layer 
and Az corresponds the width of the transition between the 
layer and the weakly magnetised atmosphere. As in section [3] 
we choose z t = 0.6 and Az = 1/60. In the following, the re- 
gion < z < z t is called the atmosphere whereas the region 
zt < z < 1 is called the layer. The x-component of the mag- 
netic field is derived from 



"■-jK-^)] 



(ll) 



The profile of horizontal field can be seen on figure [2ja), for 
e = 0.15. 

We define the pitch angle as being a = a.tMi(B y /B x ) so 
that a varies between and 7r/2, depending on the value of e. 
In the present setup, the horizontal field is six times weaker in 
the atmosphere than in the layer. The actual amplitude of the 
field in the atmosphere is rather unimportant and we conducted 
additional simulations varying this parameter. The results were 
qualitatively similar apart from the cases where the vertical gra- 
dients of the x component of the field become strong enough 
to undergo interchange instability, which makes the analysis of 
the results difficult. We therefore choose to focus on the simpler 
case where B x is large enough to have a dynamical effect on 



f http://www.vapor.ucar.edu/ 
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Figure 3. Growth rate s of the instability versus horizontal wave number 
k x for the initial configuration presented on figure [2] The results are 
presented for different angles a from (unidirectional field) to tt/2 (the 
field lines are perpendicular). We stress the case a = 7r/4 as it is the 
largest value considered numerically in the following. 



the development of the instability but small enough in order to 
remain stable to interchange modes during the simulations. We 
would like to stress that this setup is different from the case pre- 
sented in section [3] where the twist was concentrated inside the 
unstable layer. Here, we consider the usual interchange instabil- 
ity altered by the presence of an oblique field in the atmosphere 
above. 



4.1 Linear analysis 

From the linear stability analysis we find that the most readily 
destabilised modes are interchange modes (k y — 0). The in- 
stability is localised at the interface between the magnetic layer 
and the magnetised atmosphere. Figure [3] shows how the growth 
rate of the instability varies as a function of the horizontal wave 
number k x for different pitch angles a. As the pitch angle is in- 
creased the growth rate of the most unstable mode decreases due 
to the increase in magnetic tension at the interface. This lower 
growth rate means that the instability will take longer to develop 
for larger values of the twist parameter a. The mode of maxi- 
mum growth also decreases as the twist is increased, however 
the wavelength of the most unstable mode remains smaller than 
the depth of the domain for all values of the twist parameter a, 
this will be discussed further in the following sections. Although 
our configuration differs, our conclusion is qualitatively similar 
to the one drawn by |Cattaneo et al.| ( |1990a| > and |Kusano et al.| 
(1998 ), i.e. as the intensity of the transverse field is increased, 
both the growth rate of the instability and the transverse wave 
number decrease. 



4.2 Non-linear 3D numerical simulations 

In this section, we discuss the result from non-linear simulations 
in 3D. The parameter that we vary here is the initial pitch an- 
gle of the field in the atmosphere contained between z = and 
z — 0.6. We first run a reference simulation where B x is initially 
zero everywhere. This corresponds to the classical interchange- 
type instability already considered in many papers (see, for ex- 
ample, |MatffiewTet^([T995j and |Wissink et al.| ( [2000] )), apart 
from the initial presence of a weak magnetic field in the atmo- 
sphere above the layer. Using this as a reference case, we exam- 
ine three different simulations that only differ by the parameter 
e as defined in equation (TT\ . The parameter e varies from e = 



(unidirectional field, a = 0) to e = 0.15 (a = 7r/4). The reason 
why we don't consider angles larger than 7r/4 will be explained 
later. For each of these initial configurations, we solve equations 
(T])-([6| using the numerical method described previously. 

We present in figure [4] vertical slices in the plane (x,z) 
where contours of the toroidal field B y are plotted at different 
times. Bright colours correspond to large values whereas dark 
tones correspond to low values of B y . The upper part of the 
figure correspond to the unidirectional field case (i.e. B x = 
initially). As already observed in section [3] the initial instabil- 
ity consists of alternation of upward and downward flows. As 
they increase in amplitude, a secondary instability, of Kelvin- 
Helmholtz type, develops in the shear region generated at the in- 
terface between the upward and downward motions. This tends 
to generate horizontal vorticity, which further mixes field lines 
with different amplitudes of B y . The end product of this type 
of simulations is qualitatively the same for different parameters: 
the initial coherent layer is disrupted by the vortical motion gen- 
erated by the secondary instability. The resulting toroidal field 
is diffuse and no intense coherent magnetic structures are ob- 
served. It is worth noting that the presence of a weak B y field 
in the upper part of the domain doesn't modify the nature of 
the instability, as can be seen comparing the top middle panel 
of figure [4] with the top middle panel of figure [T] Note however 
that, as noted by Vasil & Brummell (2009 ), one way to damp 
the secondary instability is to consider very large values of the 
magnetic Prandtl number. By doing so, the secondary instabil- 
ity is damped by strong viscosity whereas the magnetic field re- 
mains coherent. It is however difficult to justify such a parameter 
regime as the magnetic Prandtl number based on molecular dif- 
fusivities is expected to be very small in the tachocline ( Gough 
12007) . 

Now consider a second simulation where the x component 
of the magnetic field is initially non zero in the atmosphere 
above the unstable layer. We first consider the case where the 
pitch angle is 7r/4 (as depicted in figure |2ja)) in order to illus- 
trate the dynamical effect of the transverse field B x . Initially, as 
seen on the lower left panel of figure [4] the nature of the insta- 
bility is unchanged, apart from a larger characteristic horizontal 
length scale, as predicted by the previous linear analysis of sec- 
tion |4.1| From the simulation, the most unstable wave number 
in the x direction is roughly twice as small in the case a = tt/4 
than in the reference case a = (see also the magnetic energy 
spectra in figure [9] below). At later times (t = 16.6, lower mid- 
dle panel of figure [4|, it is apparent that the non-linear develop- 
ment of the instability is also modified due the magnetic tension 
existing at the interface between the unstable layer and the atmo- 
sphere above. Figure [4] shows an alternation of strong and weak 
down flows due to the merging of small magnetic structures into 
larger ones. The reduction of the growth rate as predicted by 
the linear analysis of section [47T| is also apparent, as the mag- 
netic fluctuations are much more spread across the layer in the 
unidirectional case. At the final time t = 37 A, the difference be- 
tween the two cases is very clear: the presence of the x field in 
the atmosphere tends to generate large coherent magnetic struc- 
tures in the x direction containing fluctuations on smaller scales 
generated during the initial stages of the instability. The toroidal 
field B y is also much less diffuse than in the unidirectional case, 
and a clear interface between field and field-free regions is ob- 
served. At later times, the magnetic and velocity fluctuations are 
too close to the upper and lower boundaries, and we stop the 
simulations to avoid spurious confinement effects. 

It is informative to look in details at one of the merging 
event, where two small magnetic structures create a larger one 
in order to overcome the magnetic tension existing at the inter- 
face. A time series of a subpart of the (x,z) plane is presented 
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t = 8.3 



t = 16.6 




Figure 4. Evolution of the instability in the (x,z) plane perpendicular to the initial field lines in the layer. On the upper row, the magnetic field lines 
are unidirectional and directed perpendicularly to the plane of visualisation. On the lower row, the field lines for z < 0.6 are orientated at 45 degrees 
with respect to the field lines inside the layer. Bright colours corresponds to strong B y field whereas weaker field is dark. Time is increasing from left 
to right: t = 8.3, t = 16.6 and t = 37 A. 






Figure 5. Details of the merging mechanism. The contours show the 
current density j 2 in a small part of the (x, z) plane for different times. 
Dark tones correspond to weak current whereas light colours correspond 
to strong current. The same colorscale is used for the different times. The 
velocity field in the plane (x, z) is also shown as arrows whose length 
depends on the velocity amplitude. 



on figure [5] We plot both the current density jf 2 = ( V x B) 2 as 
contours along with arrows representing the velocity field in this 
particular vertical plane. Initially, a strong down flow in the mid- 
dle of the panel is detected, dragging down the magnetic field 
By . This vertical down flow must also advect the weak B x field 
existing in the atmosphere, therefore working against magnetic 
tension. At time t = 20.8 (third panel), it is apparent that the 
downward flow is now weaker and disconnected from the upper 
atmosphere. The secondary instability still occurs as a vortex 
has been created inside the layer. However, the vorticity is not 
continuously fed by the down flow. The final stage consists of a 
unique magnetic structure containing small scale perturbations 



slowly diffusing away. This reconnection process is clearly visi- 
ble in the current density, where two rising structures are pushed 
together due to magnetic tension, leading to the formation of one 
large current sheet at the interface. 

Finally, and before comparing quantitatively the different 
simulations, we present in figure [6] a volume rendering of the 
magnetic energy B 2 /2. Again, we compare the unidirectional 
case a = with the case a = 7r/4, at the same time t = 24. 
Without any initial transverse field B x , the results are similar 
to what was already observed by Matth ews et al.| ( |1995| ) and 
Wissink et al. (2000). The initial linear instability is purely two- 
dimensional and corresponds to the interchange of horizontal 
straight field lines. Eventually the interaction of counter-rotating 
vortices generated by the secondary instability leads to the arch- 
ing of the magnetic structures. Note that, although these mag- 
netic structures are clearly observed on figure [6] they are very 
weak and lack twist. The use of isosurfaces, as in Wissink et al. 
( 2000), can be misleading as it gives the impression the magnetic 
structures are coherent with a sharp interface, which is not the 
case. We therefore expect that the presence of overshooting con- 
vective motions at the interface between the convective and ra- 
diative zones will inevitably disrupt these small-scale magnetic 
features. In the a = tv/4 case however, we observe larger mag- 
netic structures, containing a significant amount of twist. Note 
that these larger structures are far from being well-defined lo- 
calised flux tubes ( Cattaneo et al. 2006). We however stress that 
the fact that our polytropic atmosphere is stably- stratified and 
the magnetic tension resulting from the presence of the oblique 
B x field tends to act against the rise of these structures. These 
limiting effects are the results of our over-simplified model, but 
more refined approaches could eventually lead to the formation 
of localised flux tubes from the break-up of a magnetic layer. 
The last panel on the right in figure[6]shows a close-up with sev- 
eral magnetic field lines in blue. As the strong B y field rises, 
twist accumulates at the apex of the magnetic structure, damp- 
ing its disruption by vortical motions and increasing its longevity 
against eventual convective motions. This is reminiscent of the 
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Figure 6. Volume rendering of the magnetic energy. The case a = is presented on the left whereas the case a = 7r/4 is on the middle. The 
visualisations are realised at the same time t = 24. The figure on the right presents a subvolume of the case a = 7r/4 where magnetic field lines have 
been added in order to enhance the twisted nature of the field. 
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scenario suggested by Choudhuri (2003 ), where the twist is gen- 
erated by the rise of strong toroidal fields inside the convective 
zone containing a poloidal field. Note however that Choudhuri] 
P003| ) considered already formed flux tubes and not a uniform 
layer of magnetic field. 

We now discuss more quantitatively the effect of the trans- 
verse field B x on the development of the instability. Figure [7] 
shows the evolution with time of several global quantities. In 
the following, the brackets (.) denote a spatial average over 
all coordinates. In particular, we define the total kinetic energy 
(pu 2 ) /2 and the total magnetic energy ( B 2 ) /2. As previously, 
we compare the results from simulations without initial B x field 
(a = 0) and with various pitch angles. After the decay of the 
initial perturbation, the kinetic energy grows exponentially for a 
few sound crossing times, as seen on the left panel of figure [7] 
Note that, as predicted by the linear analysis of section |4~T) the 
growth rate decreases when a increases. 

After approximately ten sound crossing times, the instabil- 
ity saturates. The evolution after this stage strongly depends on 
a. For the case a = 0, the kinetic energy decreases with time. 
In this case, most of the initial potential energy contained in the 
initial condition is released during the linear phase, leading to 
a decaying non-linear regime. For a — 7r/4, much less kinetic 
energy is initially released due to the magnetic tension at the in- 
terface, leading to a growth of the kinetic energy with time up 
to the end of the simulation. Of course, as no external source of 
energy is included in the model, we expect the kinetic energy to 
ultimately decay in all cases. One of the important steps of the 



magnetic buoyancy instability is the generation of strong hor- 
izontal vortices by the secondary Kelvin-Helmholtz instability. 
It is clear that the x field, if strong enough, will tend to damp 
this generation of vorticity. In order to create strong down and 
up flows, it is now necessary to work against the tension of the 
weak magnetic field B x . The evolution of the total enstrophy 
(u? 2 ) in the numerical domain is plotted versus time on figureP7j 
where u? = V x u. In all cases, we observe a rapid increase of 
the enstrophy with time. As the instability saturates, the growth 
rate of the enstrophy decreases. As for the kinetic energy, we 
expect the enstrophy to ultimately decay. The main effect of the 
transverse field is to reduce the production of enstrophy by the 
secondary instability, so that the enstrophy is always smaller in 
the cases where a / than in the case where a — 0. This leads 
to a reduction of the viscous dissipation, therefore reducing the 
rate at which the kinetic energy is dissipated. 

The two others panels in figure [7] show the evolution with 
time of the total magnetic energy and current density (j 2 ), 
where j = V x B. In all cases, the magnetic energy decays 
whereas the current density increases with time. As a increases, 
the magnetic energy decreases less rapidly. This is due to the 
damping of the vortical motions that reduces the mixing of the 
magnetic field lines and the creation of dissipative currents. Note 
that the current is initially slightly larger when a / 0, as ex- 
pected from the initial conditions. 

Following Stone & Gardiner (2007b ), it is useful to quan- 
tify the amount of mixing due to the instability. Instead of defin- 
ing the mixing in terms of density, as for the Rayleigh-Taylor 
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Figure 8. Mixing parameter S(z) as defined by equation fT3) versus 
depth at time t = 30. The thin black line corresponds to the initial 
profile at t = 0. 



instability, we define the fraction of strong horizontal field as 

Bh — Bh{z — 0) 



fs 



B H (z = l)-B H (z = 0) 



(12) 



where Bh — \JB% + B% . Bh{z — 0) corresponds to the ini- 
tial amplitude of the horizontal field at the upper boundary, and 
can be considered as the weak field. Bh(z = 1) is the initial 
horizontal field at the lower boundary and is therefore the strong 
field. The fraction of weak field is simply f w = 1 — f s . The 
amount of mixing in the system can be estimated by the follow- 
ing quantity 



e(z) = 4(f a f w ) 



(13) 



which is equal to when the averaged horizontal field is equal 
to its initial value at the vertical boundaries, and 1 if the field is 
perfectly mixed and is equal to the average between these two 
values. 

Figure [8] shows the mixing parameter G as a function of 
z for the four cases considered at t = 30. The thin black line 
corresponds to the value of G derived from the initial condition, 
which is nearly the same in all four cases. As time increases, the 
two fronts of mixing propagate vertically. It is clear that due to 
the magnetic tension at the interface, the overall mixing of the 
upper front is reduced. The vertical extent of the mixing region 
is reduced as a increases. This result can be explained by the 
reduced growth rate on the initial linear instability (see figure 
[3} in the presence of a transverse field and by the damping of 
the secondary instability leading to weaker vortical motions (see 
figure [7]). In the unidirectional case, the potential energy is re- 
leased to counter the stable stratification, leading to the rise of 
the toroidal magnetic field. When a / 0, the magnetic tension 
also acts against the instability. This is a clear limitation of our 
simple model, and will be further discussed in section [5] Note 
finally that figure [8] also shows that the lower front is roughly 
unchanged by the presence of transverse field. 

The generation of large-scale coherent magnetic structures 
is clearly visible in the magnetic energy spectra defined by 



E M (k x ) = ^^2^B(k x ,y,z) • B*(k x ,y,z) 



(14) 



where B(k x ,y, z) is the one-dimensional Fourier transform in 
the x direction of B(x, y, z), whereas the star denotes the com- 
plex conjugation. We focus here on the magnetic energy spectra 
in the direction perpendicular to the initial strong field (the re- 
sults in the y direction are rather similar between the different 
cases). 
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Figure 9. Magnetic energy spectrum EM(k x ) as defined by equation 
( fit) at four different times t = 2 (light tones), t = 12.5, t = 24 and 
t = 32 (dark tones), (a) a = and (b) a = tt/4. The results are 
obtained by spatially averaging over y and z. 



Figure [9] shows the magnetic energy spectrum for a = 
and a — tt/4, at four different times spanning from the initial 
linear instability (t = 2) to the late stage of the nonlinear phase 
(t = 32). Initially, the most unstable horizontal wave number 
corresponds to the linear results already discussed in section pTT] 
The most unstable wave number is roughly k x /27v w 14 for 
a = and k x /2jr w 8 for a = 7r/4, which can be compared 
to the results presented in figure [3] The time evolution of the 
magnetic energy spectra can be described in two phases. First, 
the small-scale instability grows exponentially from t « 2 to 
t ~ 12. The effect of the transverse field is to reduce the initial 
linear growth rate, as already discussed in figures [3] and|7] From 
t « 12, the evolution of the instability is constrained by the mix- 
ing due to the vortical motions. Here, the effect of the transverse 
field is to reduce the direct cascade of magnetic energy toward 
small scales, as we observe steeper energy spectra for a = tt/4 
than for a = 0, especially at early times. Finally, during the late 
stage of the instability, we observe a decay of the magnetic en- 
ergy at small scales whereas the magnetic energy at large scales 
remains constant for a = 0. In contrast, the magnetic energy is 
still growing at large scales for a = tt/4. The effect of the trans- 
verse field is indeed clearly visible on the small wave numbers 
region k x /27v < 3. For a = tt/4, a clear peak is observed at 
large scales, which corresponds to the coherent magnetic struc- 
tures already observed on figures |4]and[6] This is also consistent 
with the larger values of the global magnetic energy observed 
on figure [7] when a/0. Without transverse field, the magnetic 
energy peaks at larger horizontal wave numbers. Note that this 
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Figure 10. Relative current helicity, Hb, as defined by equation fl5) 
versus time. 



difference cannot be attributed to linear effects, as the most un- 
stable wave numbers as predicted by the linear analysis are much 
larger (see figure[3]). During the non-linear phase, the small-scale 
features (i.e. k x /2it > 4) are roughly independent of a, and we 
observe a continuity of magnetic scales up to the resistive scale. 
There is no clear scaling in our case since the diffusivities do not 
allow for a fully-developed turbulent state. 

Finally, the amount of twist in the magnetic field lines, as 
observed in figure [6] can be estimated looking at the current he- 
licity, defined by (j • B). Since both magnetic energy and cur- 
rent density are time-dependent (see figure [7}, we consider the 
relative current helicity defined by 



He 



(3-B) 



(15) 



Note that given our initial configuration, the amount of current 
helicity depends on a. We show on figure^)] the evolution with 
time of H.b for different values of a. The effect of the initial 
condition is clearly apparent, as 1~Lb = when a — 0, whereas 
1~Lb initially increases with a. Without initial transverse field, 
the current helicity remains statistically negligible, although we 
observe a very weak positive correlation, which might be due 
to the interaction between the vorticity and the rising magnetic 
field. However, as soon as a / 0, we observe an increase of the 
relative current helicity from its small initial value. This is partic- 
ularly visible in the a = ty/4 case. After a transient stage where 
1-Lb increases and then decreases, we observe a monotonic in- 
crease of the relative current helicity up to the end of the simula- 
tion. Note that Tl b continues to increase even after the saturation 
of the initial linear phase, which roughly corresponds to t ~ 20 
(see figure [7}. Again, since the atmosphere is stably stratified, 
we expect the magnetic field to ultimately stop rising and the 
current helicity to ultimately decay. If the entropy gradient was 
nearly negligible, as expected in the convective zone, the mag- 
netic structures observed in our simple model could continue to 
rise unimpeded, eventually becoming kink-unstable ( |Fan et ah] 
1999). The increase in the relative current helicity observed here 
is an indication that the twist observed in active regions could be 
due to the interaction between a rising toroidal structure and the 
weaker field inside the convective zone, as suggested by |Choud^| 
|huri| p003 ). Of course, due to numerical constraints, the large 
diffusivities considered in this paper do not allow for the forma- 
tion of isolated structures. 



5 Discussion 

In this work, we sought to examine the instability of a layer of 
horizontal magnetic field in a polytropic atmosphere, where the 
direction of the field lines is depth-dependent. In section [5] we 
examined the instability of a layer of horizontal magnetic field 
in a polytropic atmosphere, where the direction of the field lines 
in the layer varied with depth. We showed that the initial idea of 
building large-scale coherent magnetic structures from the buoy- 
ancy instability of a sheared magnetic layer, as studied by |Cat-| 
|taneo et aT] ( |1990aj ) and Kusano et al. (1998), seemed limited to 
two-dimensional geometry. The same configuration drastically 
changed in three dimensions, where the interchange instability 
was able to occur in an oblique direction. We stress that it could 
be possible to find a configuration of this type that strongly mod- 
ifies the nature of the magnetic buoyancy instability, however 
this was not the aim of the paper. 

In section [4] we introduced another initial configuration, 
where the twist was concentrated at the interface between a 
strong layer of horizontal magnetic field and a weakly mag- 
netised atmosphere above. The weakly magnetised atmosphere 
was justified by the presence of convection, which is expected 
to contain non-negligible magnetic fluctuations, both due to the 
redistribution of large-scale magnetic fields, and due to local 
small-scale dynamo action. The presence of a strong unidirec- 
tional field was justified by the presence of radial differential ro- 
tation in the tachocline. We presented a linear stability analysis 
and numerical simulations of the magnetic buoyancy instability 
in sections 14.11 and 14.21 We have shown that the initial inter- 
change instability is only weakly modified by the presence of 
the transverse magnetic field. During the non-linear phase of the 
instability, the small-scale magnetic structures merged together 
to overcome magnetic tension, leading to the formation of mag- 
netic structures on scales larger than those of the initial linear in- 
stability. These structures were shown to be more coherent than 
those in the unidirectional case, since the production of vorticity 
by the secondary instability was reduced. In addition we found 
a significant amount of twist was generated due to the rise of the 
toroidal field in the weakly magnetised atmosphere. This mech- 
anism could provide the initial twist necessary for the magnetic 
structures to rise coherently through the convective zone. 

It is however important to stress the limitations of the 
present model. In order to limit the number of parameters, we 
have considered two unidirectional fields with a current sheet at 
the interface. One could argue that the magnetic field in the con- 
vective zone is interacting with small-scale motions and is cer- 
tainly not unidirectional as assumed here. Although it remains 
to be verified using a more refined model, we expect that locally 
the toroidal field will behave similarly in the presence of small- 
scale magnetic fluctuations. 

There are two stabilising effects in the present model, the 
stably stratified atmosphere (which is necessary to suppress con- 
vective motions) and the magnetic tension generated by the 
twisted atmosphere. The tension becomes important as the at- 
mosphere is compressed by the rising structures, and the toroidal 
magnetic field will eventually stop rising. We have reached this 
regime in some of our simulations where the pitch angle is equal 
or larger that 7r/4. The stabilising effect of the tension is clearly 
overestimated by the present model, as the tension exists every- 
where at the interface between the rising layer and the atmo- 
sphere. In a more realistic situation, the convective zone would 
be filled with small- scales magnetic structures. Locally these 
structures could have a similar effect to the global tension con- 
sidered in our simulations, although we didn't check it numeri- 
cally here. We plan in a forthcoming paper to look at a config- 
uration where the global magnetic tension is removed at some 
depth, leading the rise and formation of local coherent twisted 
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magnetic structures. Another possibility would be to extend the 
work of |Fan et al.1(2003) and |Jouve & Brun| ( |2009] ), assuming 
that the convective zone is filled with magnetic perturbations. 
This would confirm whether the generation of twist, via the in- 
teraction of a rising flux tube and small-scale magnetic perturba- 
tions, is sufficient for the flux tube to remain coherent as it rises 
through the convective zone. Such simulations are underway in 
spherical geometry using the ASH code (Pinto & Brun 2012 ). 

We conclude by arguing that the merging mechanism il- 
lustrated by this simple model could have implications on the 
formation of large-scale magnetic structures from a deep-seated 
toroidal field. The twist observed in active regions could there- 
fore be attributed to a progressive accumulation of local twist as 
the toroidal field rises through the magnetised convective zone. 
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